// @HEADER
//*********************************************************************//
//  SiYuan: A numerical PDE solver                                     //
//  Copyright (2022) YUAN Xi                                           //
//  This Software is released under the BSD 2-Clause license detailed  //
//  in the file "LICENSE" in the top-level SiYuan directory            //
//*********************************************************************//
// @HEADER

#ifndef _IMPEDANCE_IMPL_HPP
#define _IMPEDANCE_IMPL_HPP

#include "Panzer_Workset.hpp"
#include "Panzer_Workset_Utilities.hpp"
#include "Panzer_IntegrationRule.hpp"

namespace SiYuan {

//**********************************************************************
// AcousticImpedance
//**********************************************************************
template<typename EvalT, typename Traits>
AcousticImpedance<EvalT, Traits> :: AcousticImpedance( const Teuchos::ParameterList& p) :
  omega_( p.get< double >("Angular Frequency") ) ,
  density_( p.get< double >("Density") )
{
  auto ir = p.get< Teuchos::RCP<panzer::IntegrationRule> >("IR");
  auto Basis = p.get< Teuchos::RCP<panzer::BasisIRLayout> >("Basis");

  auto zn = p.get< Teuchos::Array<double> >("Impedance");
  z_ = std::complex<double>(zn[0], zn[1]);

  space_dim = Basis->dimension();
  num_cub_points = Basis->numPoints();
  num_basis = Basis->cardinality();
  num_cell = Basis->numCells();
  basis_name = Basis->name();

  tempz = omega_*density_ / z_;
  
  std::string dof_name = p.get<std::string>("DOF Name");
  r_pressure_ = PHX::MDField<const ScalarT>(dof_name, Basis->functional);
  this->addDependentField( r_pressure_ );
  dof_name.replace(dof_name.end()-5, dof_name.end(), "_Img");
  i_pressure_ = PHX::MDField<const ScalarT>(dof_name, Basis->functional);
  this->addDependentField( i_pressure_ );

  std::string flux_name = p.get< std::string >("Name");
  Teuchos::RCP<PHX::DataLayout> data_layout = ir->dl_scalar;
  std::string rflux= flux_name+ "_Real";
  r_flux_ = PHX::MDField<ScalarT>(rflux, data_layout);
  this->addEvaluatedField( r_flux_ );
  std::string iflux= flux_name+ "_Img";
  i_flux_ = PHX::MDField<ScalarT>(iflux, data_layout);
  this->addEvaluatedField( i_flux_ );

  this->setName(flux_name);
}

//**********************************************************************
template<typename EvalT, typename Traits>
void AcousticImpedance<EvalT, Traits> :: postRegistrationSetup(
  typename Traits::SetupData  sd,
  PHX::FieldManager<Traits>&  fm)
{
  this->utils.setFieldData(r_pressure_,fm);
  this->utils.setFieldData(i_pressure_,fm);
  basis_index = panzer::getBasisIndex(basis_name, (*sd.worksets_)[0], this->wda);
}

//**********************************************************************
template<typename EvalT, typename Traits>
void AcousticImpedance<EvalT, Traits> :: evaluateFields(
  typename Traits::EvalData  workset)
{ //std::cout << "RES:" << PHX::print<EvalT>()  << std::endl;
  const Teuchos::RCP<const panzer::BasisValues2<double> > bv = this->wda(workset).bases[basis_index];
  for (std::size_t cell = 0; cell < workset.num_cells; ++cell) {
    for (std::size_t qp = 0; qp < num_cub_points; ++qp) {
      ScalarT Pr=0.0;
      ScalarT Pi=0.0;
      for (std::size_t basis = 0; basis < num_basis; ++basis) {
        Pr += r_pressure_(cell,basis)*bv->basis_scalar(cell,basis,qp);
        Pi += i_pressure_(cell,basis)*bv->basis_scalar(cell,basis,qp);
      }
      r_flux_(cell,qp) = - tempz.real()*Pi - tempz.imag()*Pr;
      i_flux_(cell,qp) = tempz.real()*Pr - tempz.imag()*Pi;
    //  std::cout << cell << ","  << qp << "," << r_flux_(cell,qp) << "," << i_flux_(cell,qp)  << std::endl;
    }
  }
}

//**********************************************************************
// AcousticVelocity
//**********************************************************************
template<typename EvalT, typename Traits>
AcousticVelocity<EvalT, Traits> :: AcousticVelocity( const Teuchos::ParameterList& p) :
  omega_( p.get< double >("Angular Frequency") ) ,
  density_( p.get< double >("Density") )
{
  auto ir = p.get< Teuchos::RCP<panzer::IntegrationRule> >("IR");
  auto Basis = p.get< Teuchos::RCP<panzer::BasisIRLayout> >("Basis");

  auto vn = p.get< Teuchos::Array<double> >("Normal Speed");
  v_ = std::complex<double>(vn[0], vn[1]);

  space_dim = Basis->dimension();
  num_cub_points = Basis->numPoints();
  num_basis = Basis->cardinality();
  num_cell = Basis->numCells();
  basis_name = Basis->name();
  
  wp = omega_ * density_;
//  std::string dof_name = p.get<std::string>("DOF Name");

  std::string flux_name = p.get< std::string >("Name");
  std::string rflux= flux_name+ "_Real";
  Teuchos::RCP<PHX::DataLayout> data_layout = ir->dl_scalar;
  r_flux_ = PHX::MDField<ScalarT>(rflux, data_layout);
  this->addEvaluatedField( r_flux_ );
  
  std::string iflux= flux_name+ "_Img";
  i_flux_ = PHX::MDField<ScalarT>(iflux, data_layout);
  this->addEvaluatedField( i_flux_ );

  this->setName(flux_name);
}

//**********************************************************************
template<typename EvalT, typename Traits>
void AcousticVelocity<EvalT, Traits> :: postRegistrationSetup(
  typename Traits::SetupData  sd,
  PHX::FieldManager<Traits>&  fm)
{
  basis_index = panzer::getBasisIndex(basis_name, (*sd.worksets_)[0], this->wda);
}

//**********************************************************************
template<typename EvalT, typename Traits>
void AcousticVelocity<EvalT, Traits> :: evaluateFields(
  typename Traits::EvalData  workset)
{ //std::cout << "RES:" << PHX::print<EvalT>()  << std::endl;
  const Teuchos::RCP<const panzer::BasisValues2<double> > bv = this->wda(workset).bases[basis_index];
  for (std::size_t cell = 0; cell < workset.num_cells; ++cell) {
    for (std::size_t qp = 0; qp < num_cub_points; ++qp) {
      r_flux_(cell,qp) = -wp * v_.imag();
      i_flux_(cell,qp) = wp * v_.real();
    //  std::cout << cell << ","  << qp << "," << r_flux_(cell,qp) << "," << i_flux_(cell,qp)  << std::endl;
    }
  }
}

}

#endif
